function Gycall(a)

global DAE

if ~a.n, return, end

md = DAE.y(a.md);

V1 = a.u.*DAE.y(a.v1);
V2 = a.u.*DAE.y(a.v2);  
V12 = V1.*V2;
V22 = 2.*md.*V2.*V2;
y = admittance(a);
g = real(y);
b = imag(y);
m2 = md.*md;

[s12,c12] = angles(a);

k1 = (c12.*g + s12.*b)./md;
k2 = (c12.*g - s12.*b)./md;
k3 = (s12.*g - c12.*b)./md;
k4 = (s12.*g + c12.*b)./md;

V1m = -2*V1.*V1./m2./md;

DAE.Gy = DAE.Gy ...
  - sparse(a.md, a.md, 1, DAE.m, DAE.m) ...
  - sparse(a.bus1, a.md, -V12.*k1./md - V1m.*g, DAE.m, DAE.m) ...
  - sparse(a.v1,   a.md, -V12.*k3./md + V1m.*b, DAE.m, DAE.m) ...
  - sparse(a.bus2, a.md, -V12.*k2./md, DAE.m, DAE.m) ...
  + sparse(a.v2,   a.md, -V12.*k4./md, DAE.m, DAE.m);
 
DAE.Gy = DAE.Gy + ...
          sparse(a.bus1, a.v1, 2.*V1.*g./m2-V2.*k1, DAE.m, DAE.m) + ...
          sparse(a.bus1, a.v2, -V1.*k1, DAE.m, DAE.m) + ...
          sparse(a.bus2, a.v1, -V2.*k2, DAE.m, DAE.m) + ...
          sparse(a.bus2, a.v2, 2*V2.*g-V1.*k2, DAE.m, DAE.m);

DAE.Gy = DAE.Gy + ...
          sparse(a.bus1, a.bus1, V12.*k3, DAE.m, DAE.m) - ...
          sparse(a.bus1, a.bus2, V12.*k3, DAE.m, DAE.m) + ...
          sparse(a.bus2, a.bus1, V12.*k4, DAE.m, DAE.m) - ...
          sparse(a.bus2, a.bus2, V12.*k4, DAE.m, DAE.m);

DAE.Gy = DAE.Gy + ...
          sparse(a.v1, a.v1, -2.*V1.*b./m2-V2.*k3, DAE.m, DAE.m) + ...
          sparse(a.v1, a.v2, -V1.*k3, DAE.m, DAE.m) + ...
          sparse(a.v2, a.v1,  V2.*k4, DAE.m, DAE.m) + ...
          sparse(a.v2, a.v2, -2*V2.*b+V1.*k4, DAE.m, DAE.m);

DAE.Gy = DAE.Gy - ...
          sparse(a.v1, a.bus1, V12.*k1, DAE.m, DAE.m) + ...
          sparse(a.v1, a.bus2, V12.*k1, DAE.m, DAE.m) + ...
          sparse(a.v2, a.bus1, V12.*k2, DAE.m, DAE.m) - ...
          sparse(a.v2, a.bus2, V12.*k2, DAE.m, DAE.m);

